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Abstract. We study the static screening in a Hubbard-hke model using fixed-node 
diffusion Monte Carlo. We find that the random phase approximation is surpris- 
ingly accurate even for metallic systems close to the Mott transition. As a specific 
application we discuss the implications of the efficient screening for the supercon- 
ductivity in the doped FuUerenes. In the Monte Carlo calculations we use trial 
functions with two Gutzwiller-type parameters. To deal with such trial functions, 
we introduce a method for efficiently optimizing the Gutzwiller parameters, both 
in variational and in fixed-node diffusion Monte Carlo. 



1 Introduction 

The random phase approximation (RPA) is widely used throughout solid 
state physics. R properly describes the screening when the kinetic energy is 
much larger than the Coulomb energy. In the strongly correlated limit, how- 
ever, it is qualitatively wrong. Little is known for the intermediate regime, 
where kinetic and Coulomb energy are comparable, and perturbative meth- 
ods fail since there is no small parameter. In such a situation quantum Monte 
Carlo is the method of choice. Our goal is to investigate the screening in a 
strongly correlated system close to the Mott transition. To be specific we 
focus on the doped FuUerenes, like KsCeo- These materials are strongly cor- 
related and close to a Mott transition [Q, but they are also superconduc- 
tors with quite large transition temperatures (Tc w 30K for RbsCgo) [@|- 
The superconductivity is driven by the coupling to intramolecular phonons. 
These phonons mediate an effective electron-electron attraction which is, 
however, counteracted by the electron-electron repulsion, which is large in 
such strongly correlated systems. Unlike conventional superconductors, in the 
doped FuUerenes this repulsion is not reduced much by retardation effects. 
Therefore efficient screening is important for reducing the electron-electron 
repulsion sufficiently to allow for an electron-phonon driven superconduc- 
tivity 1^,1). Using quantum Monte Carlo we find that even for quite strong 
correlations the RPA gives a surprisingly accurate description of the static 
screening on the metallic side of a Mott transition until the system is close 
to the transition Q. Besides the immediate consequences this result has for 
our understanding of the superconductivity in the doped FuUerides, it should 
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also have quite general implications for the physics of systems close to a Mott 
transition. 

We start by introducing the model Haniiltonian used to describe the 
doped FuUerenes. We briefly discuss the Lowdin downfolding technique that 
gives a general prescription for constructing low-energy Hamiltonians. Next 
we introduce the screening problem and discuss the results of the Monte Carlo 
calculations. Finally we review the quantum Monte Carlo methods used. We 
discuss the choice of the trial wave function and introduce a very efhcient 
method for optimizing the Gutzwiller-type parameters. 

2 Model Hamiltonian 

FuUerites are crystals made of Cgo molecules on an fee lattice. They are char- 
acterized by very weak inter-molecular interactions. Therefore the discrete 
molecular levels merely broaden into narrow, well separated bands The 
valence band originates from the lowest unoccupied molecular orbital, which 
is a 3-fold degenerate tiu orbital. In doped FuUerenes, there are alkali atoms 
sitting in the space between the Ceo molecules. They do not affect the band 
structure around the Fermi energy very much. Only the filling of the ti^ band 
changes, since each alkali atom donates its valence electron. Hence for KsCgo 
the (3-fold degenerate) tiu band is half-filled. We are mainly interested in 
the properties of these valence-band electrons. To simplify the description 
we therefore want to get rid of the other bands. They can be projected out 
by the Lowdin downfolding technique ||]. The basic idea is to partition the 
Hilbert space into a subspace that contains the degrees of freedom that we 
are interested in (in our case the 'ti^-subspace') and the rest of the Hilbert 
space: TC = TIq © Hi. We can then write the Hamiltonian of the system as 

where Ha is the projection of the Hamiltonian onto subspace Tii, while the 
Hij {i / j) contain the hybridization matrix elements between the two sub- 
spaces. Writing Green's function G = [E ~ H)~^ in the same way, we can 
calculate the projection of G onto Tig 0: 

Goo = [e- [Hqo + Hqi {E ~Hn)-^Hio]) . (2) 

We see that the physics of the full system is described by an effective Hamil- 
tonian Hcf[{E) that operates on the subspace Hq only. We have, however, 
to pay a price for this drastic simplification: the effective Hamiltonian is en- 
ergy dependent. In practice one approximates it with an energy-independent 
Hamiltonian iJoff(£'o). This works well if we are only interested in energies 
close to Eq. In solid Ceo we have the fortunate situation that the bands re- 
tain the character of the molecular orbitals, since the hybridization matrix 
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Fig. 1. Band structure (tiu band) of solid Ceo (fee) (a) as calculated ab initio using 
the local density approximation M and (b) using a tight-binding Hamiltonian with 
only tiu orbitals 



elements are small compared to the energy separations of the orbitals. In fact 
we can neglect the other bands altogether and get the hopping matrix ele- 
ments Un.jn' between the tiu orbitals n and n' on molecules i and j directly 
from a tight-binding parameterization [^|j9|. To demonstrate how well this 
works, Fig. [l] shows the comparison of the ab initio ii„ band structure with 
the band structure obtained from the tight-binding Hamiltonian with only 
tiu orbitals. 

A realistic description of the electrons in the ti„ band also has to include 
the correlation effects which come from the Coulomb repulsion of electrons in 
tiu orbitals on the same molecule. The resulting Hamiltonian which describes 
the interplay of the hopping of electrons and their Coulomb repulsion has the 
form 

H = tinjn' c\^a'^jn'a + ^ ninalT-in'a'- (3) 

{ij) nn' a i (ncr)<(n'cr') 

The on-site Coulomb interaction U can be calculated within density func- 
tional theory It is given by the increase in the energy of the tiu level per 
electron that is added to one molecule of the system. It is important to avoid 
double counting in the calculation of U . While the relaxation of the occupied 
orbitals and the polarization of neighboring molecules have to be included in 
the calculation, excitations within the tiu band must be excluded, since they 
are contained explicitly in the Hamiltonian (||). The results are consistent 
with experimental estimates | pl] , p^ : U ~ 1.2 — 1.4 eV. For comparison, the 
width of the tiu band is in the range W « 0.5 — 0.85 eV . To properly describe 
KaCgo the effect of the orientational disorder |^,|l3j of the Cgo molecules in 
the crystal are built into the hopping matrix elements tinjn' ■ Multiplet ef- 
fects are not included, since they tend to be counteracted by the Jahn- Teller 
effect, which is also neglected. 
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In K3C60 the system has three electrons per molecule. In the limit of weak 
correlations {U — 0), this corresponds to a metal with a half-filled conduction 
band. In the atomic limit ([/ — > 00) the Coulomb energy dominates, forcing 
every molecule to be occupied by exactly three electrons, and suppressing 
any hopping. This is a Mott insulator. We therefore expect a metal-insulator 
transition for some finite value of the Coulomb interaction U . For the model 
Hamiltonian (|^) with parameters describing K3C60 it occurs for U « 1.5 — 
1.75 eV [0. Given the estimates for the true value of U, K3C60 is therefore 
close to the Mott transition, in a correlated metallic state. 

3 Screening of a Point Charge 

We now investigate how efficient the screening is in a strongly correlated 
system like KsCeo- To be specific we analyze how a test charge q sitting on 
one molecule is screened by the conduction electrons in the ii„ band. To 
describe the influence of the test charge situated on the molecule at site c we 
include an additional term 



in the Hamiltonian (^. Determining the electron density at site c for the 
system without test charge and for the system with a finite q we find the 
screening 

An nc{0)-nc{q) 

= • (5) 

q q 

Let's first discuss the screening in the RPA. In the random phase approx- 
imation it only costs kinetic energy to screen the test charge. In the limit 
where a typical Coulomb integral U is large compared with the band width 
W, the kinetic energy cost of screening is relatively small compared with the 
potential energy gain. Therefore, within the random phase approximation, 
the screening is very efficient for large U. This means that as the test charge 
q is introduced, almost the same amount of electronic charge moves away 
from the site: An « q for large U (see Fig. ||). The random phase approxi- 
mation neglects, however, that when an electron leaves a site it has to find 
another site with a missing electron or there is an increase in Coulomb energy 
of the order of U. Thus the RPA is accurate for small values of U/W, while 
it is qualitatively wrong for large U/W . It is not clear what happens when 
Coulomb energy U and band width W are comparable. 

To address this question we have performed quantum Monte Carlo calcu- 
lations for the combined Hamiltonian 




(4) 
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Fig. 2. Screening charge An on the site of the test charge (g — 0.25 e) as a function 
of U/W, extrapolated to infinite cluster size. The full curve shows the screening 
charge in the RPA, obtained from Hartree calculations for the Hamiltonian The 
crosses with errorbars give the results of the QMC calculations. The RPA screening 
remains rather accurate up to ?7/VF ~ 2, but fails badly for larger values of U/W. 
The screening is very efficient for U/W ~ 0.5 — 2.0. 



Since we are interested in the linear response, we should calculate the effect 
of an infinitesimally small test charge q. Because of the statistical error in a 
Monte Carlo calculation it is, however, difficult to determine the response to 
a small perturbation. To get a good signal-to-noise ratio, we would therefore 
like to use as large a test charge as possible. To estimate how large we can 
make q and still be in the linear response regime, we have performed Lanc- 
zos calculations for a small system of 4 molecules, where exact calculations 
are possible. Checking the response for different test charges we find that for 
q < 0.25 e the response is practically linear. The quantum Monte Carlo cal- 
culations were then performed for large clusters of Nnio\ =32, 48, 64, 72, and 
108 molecules, corresponding to systems with 96, 144, 192, 216, and 324 elec- 
trons, and q = 0.25. To extrapolate the results Z\n(A^moi) to infinite systems 
size we used a finite-size scaling of the form An = Z\n (iVmoi) + ct/N^oi- The 
results are shown in Fig. ||. For rather small values of U/W (~ 0.5 — 1.0), 
the RPA somewhat underestimates the screening. Such a behavior is also 
found in the electron gas |Q. For intermediate values of U/W 1.0 — 2.0) 
the RPA gives surprisingly accurate results. For larger U/W the screening 
rapidly breaks down, and the RPA becomes qualitatively wrong, as discussed 
above. 

The efficient screening almost up to the Mott transition which occurs for 
U/W ~ 2.5 O has profound implications for the superconductivity in the 
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alkali doped FuUerenes. We will only give qualitative arguments here, a more 
detailed discussion can be found in |^]. In BCS theory the superconducting 
transition temperature is given by 



where A = N{0) V describes the electron-phonon coupling, which mediates 
the effective electron-electron attraction, /i* = N{0) U* is the Coulomb pseu- 
dopotential that describes the repulsive Coulomb interaction between elec- 
trons, and A^(0) is the density of states at the Fermi energy. Since the electrons 
couple to intramolecular phonons, the coupling constant is to a good ap- 
proximation a molecular property. Therefore increasing the density of states 
N(0) will increase the electron-phonon coupling A. This can be achieved by 
increasing the lattice constant of the solid FuUerene. Taking the Cgo molecules 
further apart will decrease the hopping matrix elements tinjn' thus narrow- 
ing the tiu band and correspondingly increasing the density of states. If also 
U* was a molecular quantity, i.e. independent of the lattice constant, and 
if V > U* which, of course, is the prerequisite for superconductivity, then 
according to ^ increasing the lattice constant would raise Tc. Such a varia- 
tion of the transition temperature with the lattice constant is indeed found 
experimentally Q . Here the variation of the lattice constant a is achieved by 
applying hydrostatic pressure (to reduce a) or by using 'larger' alkali metals 
like Rb and Cs to increase the lattice constant compared to K3C60 (chemical 
pressure). This suggests that by inserting even more bulky ions or molecules, 
like NH3, into the doped FuUerenes, could be increase Tc even further. As 
it turns out, however, Tc rapidly decreases with a when the lattice constant 
becomes too large, and eventually superconductivity disappears. This drop 
in Tc can be understood as a natural consequence of the breakdown of the 
screening close to the Mott transition. In the doped FuUerenes U* is essen- 
tially given by U {1 — An/q), where U is the unscreened Coulomb matrix 
element (cf. the interaction term in the Hamiltonian) and An/q describes 
the screening by the ii„ electrons. U is practically independent of the lattice 
constant, while the screening efficiency changes with the band width as shown 
in Fig. ||. For U/W in the region 1.0 — 2.0 the screening is almost RPA-like 
and does not vary very much, which means that in this region U* is almost 
independent of the lattice constant a. In this region Tc thus increases with 
a. For even smaller band width, or correspondingly larger lattice constants, 
the screening rapidly becomes inefficient. Now U* increases with the lattice 
constant, leading to a decrease of A — /x*. Hence for a too large, Tc rapidly de- 
creases with a and eventually vanishes. Since the screening only breaks down 
close to the metal-insulator transition, we have the interesting situation that 
Tc peaks close to a Mott transition! 

The efficient screening found in the Monte Carlo calculations also explains 
why the coupling to the alkali-phonons is weak. Each Cgo molecule is sur- 
rounded by 14 alkali ions with very weak force constants. When an electron 
arrives on a molecule one would therefore expect that the surrounding alkali 
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ions respond strongly. This was, however, not confirmed by experiment |L5| . 
That result can be naturally understood as an effect of the efficient screen- 
ing: When an electron arrives on a molecule other electrons leave, effectively 
leaving the molecule almost neutral. The alkali ions then only see a small 
change in the net charge and therefore respond weakly, leading to a small 
electron-phonon coupling. But being molecular crystals, doped FuUerenes 
also have intramolecular phonons. Some of those intramolecular phonons shift 
the tiu levels in such a way that the center of gravity of the energy levels is 
not changed. These are the modes that are not screened by the transfer of 
charged. They are therefore the modes that drive superconductivity in the 
Fullerides. 



4 Quantum Monte Carlo 

We now turn to the question of how the results shown in Fig. ^ were obtained. 
To keep the notation simple we will discuss the different methods for a simple 
Hubbard model (only one orbital per site, next neighbor hopping matrix 
elements tij = —t): 

H = -tY, +UY1 "*T"a • (8) 

The generalization to Hamiltonians like (|^) is straightforward. 

The first step in the quantum Monte Carlo approach is to identify a trial 
function fpT- For the Hubbard model that function should balance the oppos- 
ing tendencies of the hopping term and the interaction: Without interaction 
(i.e. for [7 = 0) the ground state of the Hamiltonian (^ is the Slater de- 
terminant that maximizes the kinetic energy. Without hopping {t = 0) 
the interaction is minimized. Since only doubly occupied sites, i.e. sites with 
riii = 1 and Uii = 1, contribute to the Coulomb energy, the electrons are 
distributed as uniformly as possible over the lattice to minimize the num- 
ber of double occupancies. A good compromise between these two extremes 
is to start from the non-interacting wavefunction but reduce the weight 
of configurations R with large double occupancies D{R). This leads (up to 
normalization) to the Gutzwiller wavefunction []l6| : 

lZ/j,(i?) = <p(i?), (9) 

with g G (0, 1] the Gutzwiller parameter. In configuration space the Coulomb 
term in the Hamiltonian is given by U D{R). Thus the Gutzwiller factor 
reflects the interaction term in the Hubbard Hamiltonian. In this spirit we 
can also construct trial functions for Hamiltonians with additional terms like 
the screening Hamiltonian (||) . We add a second Gutzwiller factor that reflects 
the interaction with the test charge q U nc{R): 

iPj,{R) = <?(i?) . (10) 

Changing the additional Gutzwiller factor h we can vary the occupation of 
the site with the test charge. 
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4.1 Variational Monte Carlo 

Since the Gutzwiller factor, like the interaction term, is diagonal in configu- 
ration space, we have to perform a sum over all configurations R in order to 
calculate the energy expectation value for a Gutzwiller wavefunction: 

_ {^t\h\^t) _ j:b.eur)i^hr) .... 

where we have introduced the local energy for a configuration R 

Since the number of configurations R grows exponentially with system-size, 
the summations in dll] ) can only be done explicitly for very small systems. For 
larger problems we use variational Monte Carlo |jl^ . The idea is to perform 
a random walk in the space of configurations, with transition probabilities 
p{R — > R') chosen such that the configurations -Rvmc in the random walk 
have the probability distribution function ^^{R). Then 

p EflvMC -^'oc(i?VMc) Y.rEio.{R)^UR) „ 

The transition probabilities can be determined from detailed balance 

W^{R) p{R R') = p{R' R) (14) 

which is fulfilled by the choice p{R R') = 1/N min[l, i?'|(i?')/^'l(-R)], 
with N being the maximum number of possible transitions. It is sufficient 
to consider only transitions between configurations that are connected by 
the Hamiltonian, i.e. transitions in which one electron hops to a neighboring 
site. The standard prescription is then to propose a transition R R' with 
probability 1/A^ and accept it with probability imn[l,^^{R')/^^{R)]. This 
works well for U not too large. For strongly correlated systems, however, the 
random walk will stay for long times in configurations with a small number 
of double occupancies D{R), since most of the proposed moves will increase 
D and hence be rejected with probability « 1 — g^'^^ 

There is, however, a way to integrate-out the time the walk stays in a 
given configuration. To see how, we first observe that for the local energy ( p^ 
the ratio of the wavefunctions for all transitions induced by the Hamiltonian 
have to be calculated. This in turn means that we also know all transition 
probabilities p(i? R'). We can therefore eliminate any rejection (i.e. accept 
with probability one) by proposing moves with probabilities 
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Checking detailed balance (14) we find that now we are sampling configura- 
tions -RvMC from the probability distribution function ^^{R) (1 — pstay(^))- 
To compensate for this we assign a weight w{R) = 1/(1 — Pstay{R)) to each 
configuration R. The energy expectation value is then given by 

^R^MC ^(^VMC) 

The above method is quite efficient since it ensures that in every Monte 
Carlo step a new configuration is created. Instead of staying in a configura- 
tion where is large, this configuration is weighted with the expectation 
value of the number of steps the simple Metropolis algorithm would stay 
there. This is particularly convenient for simulations of systems with strong 
correlations: Instead of having to do longer and longer runs as U is increased, 
the above method produces, for a fixed number of Monte Carlo steps, results 
with comparable error estimates. 



Correlated sampling So far we have only specified the form of the trial 
function ^t- The goal of a variational calculation is now to identify the pa- 
rameters that result in the best trial function. A criterion for a good trial 
function is e.g. a low variational energy. To find the wavefunction that mini- 
mizes the variational energy we could perform independent VMC calculations 
for a set of different trial functions. It is, however, difficult to compare the 
energies from these calculations since each VMC result comes with its own 
statistical errors. This problem can be avoided with correlated sampling JTsf . 
The idea is to use the same random walk for calculating the expectation value 
of all the different trial functions. This reduces the relative errors and hence 
makes it easier to find the minimum. 

Let us assume we have generated a random walk Rvmc using as the 
trial function. Using the same random walk, we can then estimate the energy 
expectation value (nSl) for a different trial function iPt, by introducing the 
reweighting factors ^{R)/^^{R): 

ER,^jHR)/^m 

with 

^ioc(i?) -5]i^^^ + C/i?(i?) (18) 



^t{R) ^t{R)/^t{R) 



U D{R) 



We notice that (also in i^ioc) the new trial function appears only in 
ratios with the old For trial functions (|^) that differ only in the Gutzwiller 
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factor this means that the Slater determinants cancel, leaving only powers 
{g / g)^^^'^ ■ Since D{R) is integer we can then rearrange the sums in (^) into 
polynomials in g/g. To find the optimal Gutzwiller parameter we then pick 
a reasonable g, perform a VMC run for ^T{g) during which we also estimate 
the coefficients for these polynomials. We can then calculate E{g) by simply 
evaluating the ratio of the polynomials. Since there are typically only of the 
order of some ten non- vanishing coefficients this method is very efficient. 

The idea of rewriting the sum over configurations into a polynomial can 
be easily generalized to trial functions with more correlation factors of the 
type r'^^^\ as long as the correlation function c{R) is integer- valued on the 
space of configurations. As a specific example of how the method works in 
practice, Fig, y shows the result of a correlated sampling run for the trial 
function ® . 



4.2 Fixed-node diffusion Monte Carlo 

Diffusion Monte Carlo allows us, in principle, to sample the true ground 
state of a Hamiltonian. The basic idea is to use a projection operator that 
has the lowest eigenstate as a fixed point. For a lattice problem where the 
spectrum is bounded £"„ G [Eq, i?max], the projection is given by 

|if = [1 - t{H - Eq)] |>f ; lip-^")) = \^t). (19) 

If r < 2/(i?niax — Eq) and \^t) has a non-vanishing overlap with the ground 
state, the above iteration converges to |!f'o)- There is no time-step error in- 
volved. Because of the prohibitively large dimension of the many-body Hilbert 
space, the matrix- vector product in ( [l9| ) cannot be done exactly. Instead, we 
rewrite the equation in configuration space 

^ ^"+'^) = - r(H - En)\R) {R\¥'"~^) (20) 

=:F{ri'M) 

and perform the propagation in a stochastic sense: is represented by 

an ensemble of configurations R with weights w{R). The transition matrix 
element F{R',R) is rewritten as a transition probability p{R R') times 
a normalization factor ■m{R',R). The iteration ( ^o| ) is then stochastically 
performed as follows: For each R we pick a new configuration R' with proba- 
bility p{R —> R') and multiply its weight by m{R' , R). Then the new ensemble 
of configurations R' with their respective weights represents Impor- 
tance sampling decisively improves the efficiency of this process by replacing 
F(R',R) with G{R',R) = {^t\R') F{R',R)/{R\^t), so that transitions from 
configurations where the trial function is small to configurations with large 
trial function are enhanced: 

^|i?'}('Z'T|i?')(-R'l'^^"+'^) - l^')G(i?',i?)('Z'r|i?)(i?|f*"^). (21) 

R,R' 
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Fig. 3. Correlated sampling for the parameters o and h in the generalized 
Gutzwiller wavefunction \^t) = g^h"" l^), cf. eqn. (|lO|), in variational (upper plot) 
and fixed- node diffusion Monte Carlo (lower plot). The plots show the energy as a 
function of the Gutzwiller parameters g and h, both as surfaces and contours. The 
calculations were done for an fee cluster of 64 molecules with 96-1-96 electrons (half- 
filled flu-band), an on-site Hubbard interaction U — 1.25 eV, and a test charge of 
q = l/4e. 

Now the ensemble of configurations represents the product S''-"-' . After a 
large number n of iterations the ground state energy is then given by the 
mixed estimator 

As long as the evolution operator has only non-negative matrix elements 
G{R' , R), all weights w{R) will be positive. If, however, G has negative matrix 
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elements there will be both configurations with positive and negative weight. 
Their contributions to the estimator ( p^ ) tend to cancel so that eventually 
the statistical error dominates, rendering the simulation useless. This is the 
infamous sign problem. A straightforward way to get rid of the sign problem is 
to remove the offending matrix elements from the Hamiltonian, thus defining 
a new Hamiltonian i?eff by 

if G(i?', i?) <0 , . 

{R'\H\R) else ^ > 

For each off-diagonal element {R'\H\R) that has been removed, a term is 
added to the diagonal: 

{R\H,s\R) ^ {R\H\R)+Y.'^t{R'){R'\H\R)I9t{R). (24) 

R' 

This is the fixed-node approximation for lattice Hamiltonians [Q. H^s is 
by construction free of the sign problem and variational, i.e. E^^ > Eq. 
The equality holds if Wt{R')/^t{R) = %{R')/MR) for all -R, R' for which 
G{R',R) is negative. 

Fixed-node diffusion Monte Carlo for a lattice Hamiltonian thus means 
that we choose a trial function from which we construct an effective Hamil- 
tonian and determine its ground state by diffusion Monte Carlo. Because 
of the variational property, we want to pick the ifr such that Eq^ is mini- 
mized, i.e. we want to optimize the trial function, or, equivalently, the effec- 
tive Hamiltonian. Also here we can use the concept of correlated sampling. 
For optimizing the Gutzwiller parameter g we can even exploit the idea of 
rewriting the correlated sampling sums into polynomials in g/g, that we al- 
ready have introduced in VMC. There is, however, a problem arising from the 
fact that the weight of a given configuration i?^") in iteration n is given by 
the product w(i?(")) = Yl7=i MR'-'\ R'''~^^) ■ Each individual normalization 
factor m(R' , R) can be written as a finite polynomial, but the order of the 
polynomial for w(i?'^")) increases steadily with the number of iterations. It is 
therefore not practical to try to calculate the ever increasing number of co- 
efficients for the correlated sampling function i?^"' (g) . But since we still can 
easily calculate the coefficients for the m{R\ R), we may use them to evaluate 
i5("^(g) in each iteration on a set of predefined values cji of the Gutzwiller pa- 
rameter. Figure ^ shows an example, again for the more general trial function 




Extrapolated estimator So far we have only considered estimators for the 
total energy. For determining the screening, however, we need to know the 
charge density ric on site c. It is no problem to calculate the expectation 
value nc(VMC) — (S'tI'^cI'^t) in variational Monte Carlo. Diffusion Monte 
Carlo gives, however, only the mixed estimator nc(DMC) = {^T\nc\^o). Since 
the density operator does not commute with the Hamiltonian, the mixed 
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Fig. 4. Screening charge An on the site of the test charge {q = 0.25 e) as a function 
of U /W, where U is the Coulomb interaction and W is the band width. Exact 
diagonahzation and QMC calculations have been performed for four molecules (12 
electrons). The figure shows that the QMC calculations are quite accurate over the 
whole range of U/W. 



estimator is different from expectation nc(QMC) = {'Fo\hc\'Po) . For a good 
choice of the trial function the true expectation value can be determined 
using the extrapolated estimator ?ic(QMC) « 2 7ic(DMC) — nc(VMC). 

To test the accuracy of this approach to the screening problem, which 
besides the extrapolated estimator also involves the fixed-node approxima- 
tion, we have compared the results of the quantum Monte Carlo calculations 
with the exact resuhs from exact diagonahzation for a small system of four 
molecules (12 electrons). As illustrated by Fig. |^, the QMC calculations for 
determining the screening charge are very accurate. 
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